% 2025.5.27
% 修改不同的电流峰值。


R_current = 1e-3;
L_current = 0.12e-3;
C_current = 84e-3;
U0_C_current_list = 100:100:1500;
num_U0 = length(U0_C_current_list);

R_test_object = 10e-3;
w0 = 1/sqrt(C_current*L_current);
alpha = (R_current+R_test_object)/2/L_current;
wd = sqrt(w0^2-alpha^2);

dt = 1e-8;      % 步长
t_final = 0.01; % 结束时间
N = round(t_final / dt) + 1; % 点数
t = linspace(0, t_final, N)';  % 生成时间列向量




%%
%转移支路
R_transfer = 10e-3;
L_transfer = 10e-6;
C_transfer = 100e-6;
U0_C_transfer = 12000;

t_final_transfer = 0.0001; % 结束时间
N_transfer = round(t_final_transfer / dt) + 1; % 点数
t_transfer1 = zeros(0.005/1e-8,1);
t_transfer2 = linspace(0, t_final_transfer, N_transfer)';  % 生成时间列向量
t_transfer3 = zeros(0.0049/1e-8,1);
t_transfer = [t_transfer1;t_transfer2;t_transfer3];

w0_transfer = 1/sqrt(C_transfer*L_transfer);
alpha_transfer = (R_transfer+R_test_object)/2/L_transfer;
wd_transfer = sqrt(w0_transfer^2-alpha_transfer^2);


if w0_transfer>alpha_transfer
     % 欠阻尼情况
     i_transfer = -U0_C_transfer/L_transfer/wd_transfer ...
     .*exp(-alpha_transfer.*t_transfer).*sin(wd_transfer.*t_transfer);
else
    error('电路非欠阻尼，解析表达式不适用此模型');
end
current_data_transfer = [t,i_transfer];
% figure;
% plot(t, -i_transfer, 'b', 'LineWidth',1.2);
% hold on;
% plot(t, -i_line,'r','LineWidth', 1.2);
% hold off;
% xlabel('时间 (ms)');
% ylabel('电流 (A)');
% title('RLC 电流响应（欠阻尼情况）');
% grid on;

%%
didt_storage = zeros(num_U0, 2);  % 每行为：U0 | I_peak
figure;
plot(t, -i_transfer, 'b', 'LineWidth',1.2);
hold on;
for cnt = 1:num_U0
    U0_C_current = U0_C_current_list(cnt);  % 当前电压值
    if w0>alpha
         % 欠阻尼情况
        i_line = -U0_C_current/L_current/wd.*exp(-alpha.*t).*sin(wd.*t);
        plot(t, -i_line,'r','LineWidth', 1.2);
        % 找出 <i> response 下峰值数值
        i_peak = max(abs(i_line));  % 默认取正最大，负峰值略过（常情况不存在负最大）
    else
        error('电路非欠阻尼，解析表达式不适用此模型');
    end
    current_data_line = [t,i_line];
    
    % 逐行计算差值（取绝对值）
    differences = abs(abs(current_data_line(1:1000001,2:2)) - abs(current_data_transfer(1:1000001,2:2)));
     % % 找出当前数组 d 的最小值 及 其对应的索引（第一次出现）
    [min_diff, index] = min(differences);
    
    % 查找所有与 min_diff 相同的索引（多个相等情况）
    all_min_indices = find(differences <=5);

    % 输出信息
    fprintf('总共有 %d 行拥有最小差值 %f\n', length(all_min_indices), min_diff);
    for k = 1:length(all_min_indices)
        i = all_min_indices(k);
        if i>50000
            t_zero = (i-500000)*1e-8;
            di_dt_zero = -U0_C_transfer/L_transfer/wd_transfer.* ...
            (-alpha_transfer.*exp(-alpha_transfer.*t_zero) ...
            .*sin(wd_transfer.*t_zero)+exp(-alpha_transfer.*t_zero).*wd_transfer.* ...
            cos(wd_transfer.*t_zero))/1e6;
            if(di_dt_zero<0)   % 第一个零点时didt为负值
                temp = differences(i);
                t_zero_next = (all_min_indices(k+1)-500000)*1e-8;
                di_dt_zero_next = -U0_C_transfer/L_transfer/wd_transfer.* ...
                (-alpha_transfer.*exp(-alpha_transfer.*t_zero) ...
                .*sin(wd_transfer.*t_zero)+exp(-alpha_transfer.*t_zero).*wd_transfer.* ...
                cos(wd_transfer.*t_zero))/1e6;
                if temp<differences(all_min_indices(k+1))
                    fprintf('%d kA峰值，过零di/dt为%.2f A/us.\n',round(i_peak/1000),di_dt_zero);
                    fprintf('行 %3d: A(i) = %.6f, B(i) = %.6f, 差值 = %f\n', ...
                        i, current_data_line(i,2), current_data_transfer(i,2), differences(i));
                else
                    if(di_dt_zero_next<0)
                        fprintf('%d kA峰值，过零di/dt为%.2f A/us.\n',round(i_peak/1000),di_dt_zero_next);
                        fprintf('行 %3d: A(i) = %.6f, B(i) = %.6f, 差值 = %f\n', ...
                            i, current_data_line((all_min_indices(k+1)),2), current_data_transfer((all_min_indices(k+1)),2), ...
                            differences((all_min_indices(k+1))));
                    else
                        fprintf('%d kA峰值，过零di/dt为%.2f A/us.\n',round(i_peak/1000),di_dt_zero);
                        fprintf('行 %3d: A(i) = %.6f, B(i) = %.6f, 差值 = %f\n', ...
                            i, current_data_line(i,2), current_data_transfer(i,2), differences(i));
                    end

                end
                peak_data(cnt, :) = [round(i_peak/1000), di_dt_zero];

                
                
                
                xlabel('时间 (ms)');
                ylabel('电流 (A)');
                title('RLC 电流响应（欠阻尼情况）');
                grid on;
                break;
             end
        end
    
    end
end
hold off;

% 绘制didt随电流峰值变化图
figure;
plot(peak_data(:,1), peak_data(:,2), 'bo-', 'LineWidth', 1.0, 'MarkerFaceColor', 'r');

% 设置坐标轴范围
axis([0 35 -1300 -500]);  % x: (100 ~ 1500V), y: (0 ~ 25A 假设最大为22A)


grid on;

xlabel('电流峰值 /kA','FontSize', 18);
ylabel('didt (A/us)', 'FontSize', 18);
title('didt随电流峰值变化关系');






